Spectra Story

# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(lmerTest)
library(emmeans)
library(merTools)
library(tidyquant)
library(ggeffects)
library(performance)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(ggdist)
library(pander)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflicts_prefer(lmerTest::lmer)

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    rect = element_rect(fill = "white", color = NA), 
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    legend.position = "bottom"))



# vectors for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
area_levels_all <- c("Northeast Shelf", area_levels)
area_levels_long_all <- c("Northeast Shelf", area_levels_long)

# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "All"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"



# Function to get the Predictions
# Drop effect fits that are non-significant  ###
get_preds_and_trendsignif <- function(mod_x){
  modx_preds <- as.data.frame(
    # Model predictions
    ggpredict(
      mod_x, 
      terms = ~ est_year * survey_area * season) ) %>% 
    mutate(
      survey_area = factor(group, levels = area_levels),
      season = factor(facet, levels = c("Spring", "Fall")))
  
    # Just survey area and year
    modx_emtrend <- emtrends(
      object = mod_x,
      specs =  ~ survey_area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
  
    # Preds with signif
    modx_preds %>% left_join(select(modx_emtrend, survey_area, season, non_zero))
  
}

Storycrafting - Northeast Shelf Spectra

I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills

This document will be split into two sections (the second potentially moving to its own doc).

Part 1: Summary of Trends

For the first section the focus is on understanding changes in the size structure. This will be done using length and bodymass spectra, large and small fish indices, and looking at how body mass size quantiles have moved over time.

This section should address: Has the size distribution changed over time, where has it changed, and some nuance around what aspects about the distribution have shifted.

Starting from Square Uno

I’d like for us to start by rebuilding the story of how the size structure has changed over time and regional differences.
As such, I would like for us to look at plots of the size spectrum slope for (1) the all-species length spectrum and (2) the Wigley species biomass spectrum from annual, seasonal and regional perspectives.
I am thinking this would be two panel plots (all-spp length, Wigley-spp biomass) with 5 subpanels each (NES, GOM, GB, SNE, MAB) that are set up to show lines by seasonal and annual time steps, much like the plot on the right of slide 7 in your most recent slide deck…except with NES added and an annual line added. Keeping the significant trend lines in is helpful. This is the first step, so when you have that ready, please share it with the group. - KMills

# Data used for Wigley estimates
wigley_in <- read_csv(here::here("Data/processed/wigley_species_trawl_data.csv"))
finfish_in <- read_csv(here::here("Data/processed/finfish_trawl_data.csv"))

Aside: Biomass Fraction in Each Community

If we need to make a decision about which community we want to use, here is the different biomass totals for each:

# glimpse(wigley_in)

# Load raw
# General tidying only, removal of strata outside our study area, Spring and Fall only
# (removes inshore and rarely/inconsistently sampled strata etc.)
# shrimps removed
trawl_basic <- read_csv(here::here("Data/processed/trawl_clean_all_data.csv"))

# Biomass totals
raw_totals <- trawl_basic %>% 
  distinct(est_towdate, cruise6, station, stratum, tow, comname,  biomass_kg) %>%   pull(biomass_kg) %>% 
  sum()

wig_biomass_total <- wigley_in %>% 
  distinct(est_towdate, cruise6, station, stratum, tow, comname,  biomass_kg) %>%   pull(biomass_kg) %>% 
  sum()

ffish_biomass_total <- finfish_in %>% 
  distinct(est_towdate, cruise6, station, stratum, tow, comname,  biomass_kg) %>%   pull(biomass_kg) %>% 
  sum()

# Make a table for display
tibble(
  "Community" = c("All (including crustaceans)", "All finfish", "Wigley Species Subset"),
  "Biomass Total" = c(raw_totals, ffish_biomass_total, wig_biomass_total)) %>% 
  mutate(`% of Total` = (`Biomass Total` / raw_totals)*100) %>% 
  gt()

Community	Biomass Total	% of Total
All (including crustaceans)	3781548	100.00000
All finfish	3708046	98.05629
Wigley Species Subset	3539719	93.60501

A. Size Spectra Trends

For the figures in part 1: non-zero trends over time have been highlighted in the appropriate plot panels with linear trends overlaid.

# Load the three regional datasets
# wigley_lenspectra <- read_csv(
#   here::here("Data/processed/wigley_species_length_spectra.csv"))  %>% 
#   mutate(survey_area = factor(survey_area, levels = area_levels))
wigley_bmspectra <- read_csv(
  here::here("Data/processed/wigley_species_bodymass_spectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))
finfish_lenspectra <- read_csv(
  here::here("Data/processed/finfish_length_spectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))

# Load Shelf
finfish_lenspectra_shelf <- read_csv(
  here::here("Data/processed/shelfwide_finfish_length_spectra.csv"))
wigley_bmspectra_shelf <- read_csv(
  here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))




# Join them together
wigley_bmspectra_all <- bind_rows(wigley_bmspectra, wigley_bmspectra_shelf) %>% 
  mutate(metric = "bodymass_spectra",
         community = "Wigley Subset",
         survey_area = factor(survey_area, levels = area_levels_all))
ffish_lenspectra_all <- bind_rows(finfish_lenspectra, finfish_lenspectra_shelf) %>% 
    mutate(metric = "length_spectra",
           community = "Finfish Community",
         survey_area = factor(survey_area, levels = area_levels_all))
# both metrics+communities
spectra_results <-  bind_rows(wigley_bmspectra_all, ffish_lenspectra_all) %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels_all),
    metric = fct_rev(metric))

Wigley Species - Bodymass Spectra

# Model linear trends in time
bmspectra_mod_wig_all <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = wigley_bmspectra_all)

# Get the predictions and flag whether they are significant
bmspectra_trend_predictions <- get_preds_and_trendsignif(bmspectra_mod_wig_all) %>% 
  mutate(model_id = "bodymass_spectra-Wigley Subset") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric),
         survey_area = factor(survey_area, levels = area_levels_all))





# Plot over observed data
# Contrast seasonal differences




# Make the plot
bmspectra_trend_p <- ggplot() +
   geom_ribbon(
    data = filter(bmspectra_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = wigley_bmspectra_all,
    aes(est_year, b, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = wigley_bmspectra_all,
    aes(est_year, b, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(bmspectra_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community, scales = "free") +
  labs(
    x = "Year",
    y = "Exponent of Body Mass Distribution",
    linetype = "Linetype",
    color = "Season",
    fill = "Season")
 
# Show plot
bmspectra_trend_p

# # Summary tables
# bmspectra_mod_wig_all %>% 
#   gtsummary::tbl_regression()   %>% 
#   add_global_p(keep = T) %>% 
#   bold_p(t = 0.05)  %>%
#   bold_labels() %>%
#   italicize_levels() %>% 
#   as_gt() %>% 
#   tab_header(title = "Wigley Community - Mass Spectra Trends")

Full Finfish Community - Length Spectra

# Model linear trends in time
lenspectra_mod_all <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = ffish_lenspectra_all)

# Get the predictions and flag whether they are significant
lenspectra_trend_predictions <- get_preds_and_trendsignif(lenspectra_mod_all) %>% 
  mutate(model_id = "length_spectra-Finfish Community") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric),
         survey_area = factor(survey_area, levels = area_levels_all))





# Plot over observed data
# Contrast seasonal differences




# Make the plot
lenspectra_trend_p <- ggplot() +
   geom_ribbon(
    data = filter(lenspectra_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = ffish_lenspectra_all,
    aes(est_year, b, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = ffish_lenspectra_all,
    aes(est_year, b, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(lenspectra_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community, scales = "free") +
  labs(
    x = "Year",
    y = "Exponent of Body Length Distribution",
    linetype = "Linetype",
    color = "Season",
    fill = "Season")
 
# Show the plot
lenspectra_trend_p

# # Summary tables
# lenspectra_mod_all  %>% 
#   gtsummary::tbl_regression()   %>% 
#   add_global_p(keep = T) %>% 
#   bold_p(t = 0.05)  %>%
#   bold_labels() %>%
#   italicize_levels() %>% 
#   as_gt() %>% 
#   tab_header(title = "Full Finfish Community - Length Spectra Trends")

Both Community Size Spectra

# Left side:
# All species  species length spectrum 
# color = season vs. annual
# facet_wrap(~survey_area)

# Right side:
# Wigley species biomass spectrum
# color = season vs. annual
# facet_wrap(~survey_area)


# Plot over observed data
# Contrast seasonal differences
(
  (lenspectra_trend_p + theme(strip.text.y = element_blank()) | 
  bmspectra_trend_p)
  ) + plot_layout(guides = "collect")

B. Median Length & Weight Trends

I also think it would be valuable to build out the story of size change further, so am thinking that plots like those of average (or median) length (all species) and weight (wigley species) for NES and the sub-regions, like those on p. 12 of the attached file above (the very first draft from late 2022) would be useful. - KMills

# Load the median weight/length data
wigley_size_df <- read_csv(here::here("Data/processed/wigley_species_size_summary.csv"))
finfish_size_df <- read_csv(here::here("Data/processed/finfish_species_length_summary.csv"))

# shelf-wide summary
finfish_lengths_shelf <- read_csv(here::here("Data/processed/shelfwide_finfish_species_length_summary.csv"))
wigley_sizes_shelf <- read_csv(here::here("Data/processed/shelfwide_wigley_species_size_summary.csv"))




# Join them
size_results <- bind_rows(
  list(
    "Finfish Community" = bind_rows(finfish_size_df, finfish_lengths_shelf),
    "Wigley Subset" = bind_rows(wigley_size_df, wigley_sizes_shelf)), 
  .id = "community"
) %>% 
  mutate(
    survey_area = factor(survey_area, levels = area_levels_all)
  )




# Get trends:
len_mod_ffish <- lmer(
  formula = med_len_cm ~ est_year * survey_area * season + (1|est_year),
  data = finfish_size_df)
# len_mod_wig <- lmer(
#   formula = med_len_cm ~ est_year * survey_area * season + (1|est_year),
#   data = wigley_size_df)
wt_mod_wig <- lmer(
  formula = med_wt_kg ~ est_year * survey_area * season + (1|est_year),
  data = drop_na(wigley_size_df, med_wt_kg))



# Get the predictions and flag whether they are significant
size_trend_predictions <- bind_rows(list(
  "med_len_cm-Finfish Community" = get_preds_and_trendsignif(len_mod_ffish),
  "med_wt_kg-Wigley Subset"   = get_preds_and_trendsignif(wt_mod_wig)
  ), .id = "model_id") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric))






# Left side: 
# Median length - all species
# color = season vs. annual
# facet_wrap(~survey_area)
size_long <- size_results %>% 
  pivot_longer(cols = c(med_len_cm, med_wt_kg), 
               names_to = "metric", 
               values_to = "val")  %>% 
  mutate(survey_area = fct_relevel(survey_area, area_levels_all),
         season = fct_rev(season))


# just length for scales
len_plot <-  size_long %>% 
  filter(metric == "med_len_cm",
         community == "Finfish Community") %>% 
  drop_na(val) %>% 
  ggplot() +
  geom_ribbon(
    data = filter(size_trend_predictions, metric == "med_len_cm", non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    #data = filter(size_long, metric == "med_len_cm"),
    aes(est_year, val, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, val, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(size_trend_predictions, metric == "med_len_cm", non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community,
             scales = "free") +
  theme(strip.text.y = element_blank()) +
  labs(
    x = "Year",
    y = "Length (cm)",
    color = "Season",
    linetype = "Linetype",
    fill = "Season")

# Right: 
# Median weight - wigley species
# color = season vs. annual
# facet_wrap(~survey_area)

# weight plots
wt_plot <- size_long %>% 
  filter(metric == "med_wt_kg") %>% 
  drop_na(val) %>% 
  ggplot() +
  geom_ribbon(
    data = filter(size_trend_predictions, metric == "med_wt_kg", non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    aes(est_year, val, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, val, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(size_trend_predictions, metric == "med_wt_kg", non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community,
             scales = "free")  +
  labs(
    x = "Year",
    y = "Weight (kg)",
    color = "Season",
    linetype = "Linetype",
    fill = "Season")

(len_plot | wt_plot) + plot_layout(guides = "collect") 

# (len_plot | wt_plot) + plot_layout(guides = "collect", widths = c(3,1.5)) & theme(legend.position = "bottom")

Notes:

There are a number of years when median weight surges in MAB. For these situations I think we likely could trace these down to exceptionally high catches of larger species like elasmobranchs.

We also don’t see a shelf-wide decline in size so I might need to check what Friedland did or prepare for a confrontation.

C. Small/Large Fish Trends (True Percentiles)

Data for part C. comes from the Wigly species subset to set percentiles based on individual bodymass.

The following figures frame “large fish” as individuals with bodymass greater than the 90th percentile for all years within that region. “Small fish” are considered smaller than the 10th percentile of the mass distribution.

The following table shows what those percentiles are for each region.

NOTE: Percentiles are determined using DescTools::Quantile() which uses numlen_adj to weight everything by abundance.

# # Get 95th percentile as threshold
# DescTools::Quantile(
#   wigley_in$length_cm, 
#   weights = wigley_in$numlen_adj, 
#   probs = c(0.05, 0.1, 0.5, 0.90, .95)) %>% 
#   t() %>% 
#   as_tibble()


# Do them separately
area_size_quants <- bind_rows(
    mutate(wigley_in, survey_area = "Northeast Shelf"),
    wigley_in) %>% 
  split(.$survey_area) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.05, 0.1, 0.5, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "survey_area") %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  arrange(survey_area)

# Show the table
area_size_quants %>% 
  mutate_if(is.numeric, round, 2) %>% 
  gt() %>% 
  tab_header(title = "Wigley Species - Body Mass Quantiles (kg)") %>% tab_spanner(label = "Percentile", columns = -survey_area) %>% 
  cols_label("survey_area" = "")

Wigley Species - Body Mass Quantiles (kg)
	Percentile
	5%	10%	50%	90%	95%
Northeast Shelf	0.01	0.01	0.08	0.93	1.61
GoM	0.01	0.02	0.14	1.11	1.68
GB	0.00	0.01	0.14	0.91	1.51
SNE	0.00	0.01	0.05	0.58	1.43
MAB	0.01	0.02	0.04	1.21	1.68

##### Large Fish Index  ####


# Assign large/small
wigley_in <- wigley_in %>% 
  left_join(select(area_size_quants, survey_area, `10%`, `90%`)) %>%   
  mutate(
    small = ifelse(ind_weight_kg < `10%`, T, F),
    large = ifelse(ind_weight_kg > `90%`, T, F))



# Get total abundance
total_abundbio <- bind_rows(
    mutate(wigley_in, survey_area = "Northeast Shelf"),
    wigley_in) %>% 
  group_by(survey_area, est_year, season) %>% 
  summarise(
    total_abund = sum(numlen_adj),
    total_bio = sum(sum_weight_kg),
    .groups = "drop")
  
  
# Get large fish bio
lf_abundbio <- bind_rows(
    mutate(wigley_in, survey_area = "Northeast Shelf"),
    wigley_in) %>% 
  group_by(survey_area, est_year, season) %>% 
  summarise(
    small_abund = sum(numlen_adj * small),
    small_bio = sum(sum_weight_kg * small),
    large_abund = sum(numlen_adj * large),
    large_bio = sum(sum_weight_kg * large),
    .groups = "drop")

# Join and divide
lfi <- left_join(total_abundbio, lf_abundbio) %>% 
  mutate(
    LFI = large_abund / total_abund,
    SFI = small_abund / total_abund)  %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  left_join(
    select(area_size_quants, survey_area, `10%`, `90%`) %>% 
      mutate(across(-survey_area, \(x) str_c(round(x, 2), " kg"))))

The following species are those which has representation in the large fish index:

You might think: “Why is a 1kg acadian redfish a”large fish”?“. That’s a consequence of abundance declining with size pretty dramatically.

big_fishes <- wigley_in %>% 
  filter(large) %>% 
  distinct(comname) %>% 
  arrange(comname) %>% 
  rename(`Common Name:` = comname) 


# Print it neatly
pander(big_fishes$`Common Name:`)

acadian redfish, alewife, american plaice, american shad, atlantic angel shark, atlantic cod, atlantic croaker, atlantic halibut, atlantic mackerel, atlantic sharpnose shark, atlantic sturgeon, atlantic torpedo, atlantic wolffish, barndoor skate, black sea bass, blackbelly rosefish, blueback herring, bluefish, bluntnose stingray, buckler dory, bullnose ray, butterfish, clearnose skate, cownose ray, cunner, cusk, fawn cusk-eel, fourspot flounder, goosefish, greater amberjack, haddock, little skate, longhorn sculpin, northern searobin, ocean pout, offshore hake, pollock, red hake, roughtail stingray, sand tiger, sandbar shark, scup, sea raven, silver hake, smooth butterfly ray, smooth dogfish, smooth skate, southern stingray, spanish mackerel, spiny butterfly ray, spiny dogfish, spotted hake, striped bass, summer flounder, tautog, thorny skate, weakfish, white hake, windowpane flounder, winter flounder, winter skate, witch flounder and yellowtail flounder

There’s really just so few “large” individuals that the percentiles sit quite small.

wigley_in %>% 
  filter(ind_weight_kg <= 5) %>% 
  ggplot(aes(x = ind_weight_g, y = ..density.., weight = numlen_adj)) + 
  geom_histogram(fill = "#535353", binwidth = 50 ) +
  #geom_density(color = "#EA4F12", linewidth = 1) +
  labs(title = "Overall Bodymass Distribution for Wigley Species",
       subtitle = "Weights over 5kg removed",
       x = "Mass (g)")

Shifting 10th & 90th Percentiles

A 90th percentile fish is smaller than I feel it has been in the past, so I wanted to plot those quantiles over time.

Indeed, the 90th & 95th percentile size for fishes in the GOM was nearly 2x what it is today.

# Do them separately
yearly_size_quants <- bind_rows(
    mutate(wigley_in, survey_area = "Northeast Shelf"),
    wigley_in) %>% 
  unite(col = "group_var", survey_area, est_year, sep = "X", remove = T) %>% 
  split(.$group_var) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "group_var") %>% 
  separate(group_var, into = c("survey_area", "est_year"), sep = "X") %>% 
  mutate(
    est_year = as.numeric(est_year),
    survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  arrange(survey_area)


# Reshape to make filtering easier
quants_long <- yearly_size_quants %>% 
  #select(-c(`5%`,`50%`, `95%`)) %>% 
  pivot_longer(cols = -c(survey_area, est_year), names_to = "Percentile", values_to = "biomass_kg") 


# Separate axes
q10_p <- quants_long %>%
  filter(Percentile == "10%") %>% 
  ggplot(aes(est_year, biomass_kg)) +
  geom_point(size = 0.8, alpha = 0.5) +
  geom_ma(n = 5, ma_fun = SMA, aes(linetype = "5-Year Moving Average"), color = "black") +
  scale_color_gmri() +
  theme(
    strip.text.y = element_blank(),
    plot.margin = margin(0,0,0,0)) +
  facet_grid(survey_area~Percentile) +
  labs(y = "Bodymass (kg)", x = "Year", linetype = "")

# Separate axes
q90_p <- quants_long %>%
  filter(Percentile == "90%") %>% 
  ggplot(aes(est_year, biomass_kg)) +
  geom_point(size = 0.8, alpha = 0.5) +
  geom_ma(n = 5, ma_fun = SMA, aes(linetype = "5-Year Moving Average"), color = "black") +
  scale_color_gmri(reverse = T) +
  theme(
    strip.text.y = element_blank(),
    plot.margin = margin(0,0,0,0)) +
  facet_grid(survey_area~Percentile) +
  labs(y = "", x = "Year", linetype = "")


# Separate axes
q95_p <- quants_long %>%
  filter(Percentile == "95%") %>% 
  ggplot(aes(est_year, biomass_kg)) +
  geom_point(size = 0.8, alpha = 0.5) +
  geom_ma(n = 5, ma_fun = SMA, aes(linetype = "5-Year Moving Average"), color = "black") +
  scale_color_gmri(reverse = T) +
  theme(plot.margin = margin(0,0,0,0)) +
  facet_grid(survey_area~Percentile) +
  labs(y = "", x = "Year", linetype = "")

(q10_p | q90_p | q95_p) + plot_layout(guides = "collect", axes = "collect_x") +
  plot_annotation(title = "Regional Size Percentile Changes")

10th percentile weights seem stable. 90th percentile body-weights have fluctuated over time.

# Map over years and areas?
quant_fillers <- yearly_size_quants %>% 
  unite(col = "year_area", sep = "X", est_year, survey_area, remove = F) %>% 
  split(.$year_area) %>% 
  map_dfr(function(x){
    tribble(
      ~"quant_range",   ~"x",            ~"y",       ~"xend",         ~"yend",
       "0-10%",   x[["est_year"]], 0,          x[["est_year"]], x[["10%"]],
       "10-20%",   x[["est_year"]], x[["10%"]], x[["est_year"]], x[["20%"]],
      #"0-20%",   x[["est_year"]], 0, x[["est_year"]], x[["20%"]],
      # Normal below
      "20-30%",   x[["est_year"]], x[["20%"]], x[["est_year"]], x[["30%"]],
      "30-40%",   x[["est_year"]], x[["30%"]], x[["est_year"]], x[["40%"]],
      "40-50%",   x[["est_year"]], x[["40%"]], x[["est_year"]], x[["50%"]],
      "50-60%",   x[["est_year"]], x[["50%"]], x[["est_year"]], x[["60%"]],
      "60-70%",   x[["est_year"]], x[["60%"]], x[["est_year"]], x[["70%"]],
      "70-80%",   x[["est_year"]], x[["70%"]], x[["est_year"]], x[["80%"]],
      "80-90%",   x[["est_year"]], x[["80%"]], x[["est_year"]], x[["90%"]]#,
      # Don't go to 100
      #"90-95%",  x[["est_year"]], x[["90%"]], x[["est_year"]], x[["95%"]]
      )

  }, .id = "year_area") %>% 
  separate(year_area, into = c("est_year", "survey_area"), sep = "X") %>% 
  mutate(
    quant_range = factor(quant_range, levels = c(
     "0-10%", "10-20%", 
     #"0-20%", 
     "20-30%", "30-40%", "40-50%", 
     "50-60%", "60-70%", "70-80%", "80-90%", "90-95%")),
    survey_area = factor(survey_area, levels = area_levels_all))
  


# Try that funky guy
ggplot(quants_long) +
  geom_ribbon(
    data = quant_fillers,
    aes(x = x, ymin = y*1000, ymax = yend*1000, fill = quant_range), 
    alpha = 1, color = "transparent", linewidth = 0.1) +
  geom_line(
    aes(est_year, biomass_kg*1000, group = Percentile), 
    linewidth = 0.2, alpha = 0.35, color = "white") +
  
  facet_grid(survey_area~., scales = "free") +
  guides(fill = guide_legend(nrow = 2, title.position = "top", title.hjust = 0.5)) +
  scale_fill_brewer(direction = 1) +
  scale_x_continuous(expand = expansion(add = c(0,0))) +
  scale_y_log10() +
  labs(
    y = "Individual Bodymass (g)", 
    x = "Year", 
    fill = "Body Size Quantile Range",
    title = "Changes in Body Size Quantiles")

LFI & SMI

The axes on these are: \[\frac{\sum{individuals_{ijs} > mass~threshold}} {\sum{individuals_{ijs}}} * 100\] for each: \(year_i\), \(region_j\), \(season_s\).

# Plot them as loing-term average

# Plot the SFI
sfi_p <- lfi %>% 
  ggplot() +
  geom_point(
    aes(est_year, SFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, SFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  theme(strip.text.y = element_blank()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "Percent of Total Abundance",
    title = "Small Fish Index\n(< 10th Percentile Bodymass)",
    fill = "Season",
    color = "Season",
    linetype = "")


# Plot the LFI
lfi_p <- ggplot(lfi) +
  geom_point(
    aes(est_year, LFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, LFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "",
    title = "Large Fish Index\n(> 90th Percentile Bodymass)",
    fill = "Season",
    color = "Season",
    linetype = "")



# Plot together
(sfi_p | lfi_p) + plot_layout(guides = "collect")

Large and Small Fish Abundance

This plots the abundance of individuals above/below the size thresholds as abundance.

# Plot them as loing-term average

# Plot the SFI
sfi_p <- lfi %>% 
  ggplot() +
  geom_point(
    aes(est_year, small_abund, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, small_abund, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_number()) +
  theme(strip.text.y = element_blank()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "Total Abundance",
    title = "Small Fish Index\n(< 10th Percentile Bodymass)",
    fill = "Season",
    color = "Season",
    linetype = "")


# Plot the LFI
lfi_p <- ggplot(lfi) +
  geom_point(
    aes(est_year, large_abund, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, large_abund, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA) +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_number()) +
  facet_grid(survey_area~., scales = "free") +
  labs(
    x = "Year",
    y = "",
    title = "Large Fish Index\n(> 90th Percentile Bodymass)",
    fill = "Season",
    color = "Season",
    linetype = "")




(sfi_p | lfi_p) + plot_layout(guides = "collect")

D. 1kg Discrete Size Breaks

Due to the skewed distribution in abundance with increasing size, using percentile cutoffs leads to a situation where “large” individuals are in practice quite small.

For the plots below I took the the 1-5kg range in body weights and broke that range into five equally spaced body size ranges.

# What percent is 5kg and smaller?
wigley_smol <- filter(wigley_in, ind_weight_kg <= 5)
smol_frac <- signif((sum(wigley_smol$numlen_adj, na.rm = T) / sum(wigley_in$numlen_adj, na.rm = T)), 3) *100

This still includes smol_frac% of the abundance for wigley species.

# Do everything together
# Make evenly spaced bins
overall_weight_ranges <- seq(
      floor(min(wigley_smol$ind_weight_kg, na.rm = T)),
      ceiling(max(wigley_smol$ind_weight_kg, na.rm = T)),
      length.out = 6)
# Make that a table for plotting
overall_range_tbl <-tibble(
      "size_group" = c("S", "S-M", "M", "M-L", "L"),
      "min" = overall_weight_ranges[1:5],
      "max" = overall_weight_ranges[2:6]
    ) %>% 
  mutate(size_range = str_c(min, "-",max, "kg"))


# Plot the breaks
overall_range_tbl  %>% 
  mutate(
    size_group = factor(size_group, levels = c("S", "S-M", "M", "M-L", "L"))) %>% 
  ggplot() +
  geom_segment(aes(x = min, xend = max, y = 0, yend = 0, color = size_group), linewidth = 14) +
  geom_point(aes(x = min, y = 0), shape = "|", size = 5) +
  geom_point(aes(x = max, y = 0), shape = "|", size = 5) +
  scale_y_continuous(expand = expansion(add = c(0.001, 0.001)), labels = NULL, breaks = NULL) +
  rcartocolor::scale_color_carto_d(palette = "OrYel") +
  labs(x = "Individual Body Mass (kg)", y = "", color = "Size Range", title = "Even Breaks in Body Mass",
       subtitle = ">5kg Individuals Removed")

If we plot these 1kg bins we can see that 0-1kg individuals dominate the region.

# Apply new breaks and get trends
wigley_labelled <- wigley_in %>% 
  mutate(
    size_group = case_when(
      ind_weight_kg < 1 ~ "S",
      ind_weight_kg < 2 ~ "S-M",
      ind_weight_kg < 3 ~ "M",
      ind_weight_kg < 4 ~ "M-L",
      ind_weight_kg <= 5 ~ "L",
      ind_weight_kg > 5 ~ NA_character_)) %>% 
  drop_na(size_group) %>% 
  mutate(size_group = factor(size_group, levels = c("S", "S-M", "M", "M-L", "L"))) %>% 
  left_join(overall_range_tbl)

#wigley_labelled %>% filter(size_group == "S") %>% glimpse()


# Get totals:
even_breaks_summary <- bind_rows(
    mutate(wigley_labelled, survey_area = "Northeast Shelf"),
    wigley_labelled) %>% 
  group_by(survey_area, est_year, season, size_group, size_range) %>% 
  summarise(
    total_size_abund = sum(numlen_adj, na.rm = T),
    total_size_bio = sum(ind_weight_kg),
    .groups = "drop") %>% 
  mutate(
    est_year = as.numeric(est_year),
    survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  arrange(survey_area) %>% 
  left_join(total_abundbio) %>% 
  mutate(
    size_group = factor(size_group, levels = c("S", "S-M", "M", "M-L", "L")),
    percent_abund = total_size_abund / total_abund,
    percent_bio = total_size_bio / total_bio,
    survey_area = factor(survey_area, levels = area_levels_all),
    season = factor(season, levels = c("Spring", "Fall")))



# Bar Chart
ggplot(even_breaks_summary, aes(est_year, percent_abund, fill = size_group)) +
  geom_col(position = "stack", color = "black", linewidth = 0.1) +
  rcartocolor::scale_fill_carto_d(palette = "OrYel") +
  scale_y_continuous(labels = label_percent(), limits = c(0,1), expand = expansion(add = c(0,0))) +
  facet_grid(survey_area ~ season, scales = "free") +
  labs(y = "% / Seasonal Abundance", x = "Year", fill = "Size Class")

Most Abundant Species

# What are the most abundant species in each group
top_species <- wigley_labelled %>% 
  group_by(survey_area, size_group, comname) %>% 
  summarise(total_abund = sum(numlen_adj, na.rm = T),
            total_bio  = sum(ind_weight_kg, na.rm = T)) 
# Most abundant
abund_table <- top_species %>% 
  arrange(desc(total_abund)) %>% 
  mutate(rank = row_number()) %>% 
  slice_head(n = 3) %>% 
  ungroup() %>% 
  select(-starts_with("total")) %>% 
  pivot_wider(names_from = size_group, values_from = comname) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels)) %>% 
  arrange(survey_area) %>% 
  select(survey_area, rank, S, `S-M`, M, `M-L`, L) %>% 
  group_by(survey_area) %>% 
  gt() %>% 
    tab_header(title = "Most Abundant Overall")

abund_table

Most Abundant Overall
rank	S	S-M	M	M-L	L
GoM
1	silver hake	spiny dogfish	spiny dogfish	spiny dogfish	atlantic cod
2	acadian redfish	haddock	haddock	atlantic cod	pollock
3	atlantic herring	white hake	white hake	white hake	spiny dogfish
GB
1	haddock	spiny dogfish	spiny dogfish	winter skate	winter skate
2	silver hake	haddock	winter skate	spiny dogfish	atlantic cod
3	spiny dogfish	winter skate	atlantic cod	atlantic cod	spiny dogfish
SNE
1	butterfish	spiny dogfish	spiny dogfish	spiny dogfish	spiny dogfish
2	silver hake	winter skate	winter skate	smooth dogfish	winter skate
3	spiny dogfish	ocean pout	smooth dogfish	winter skate	goosefish
MAB
1	butterfish	spiny dogfish	spiny dogfish	spiny dogfish	spiny dogfish
2	spotted hake	smooth dogfish	smooth dogfish	smooth dogfish	smooth dogfish
3	scup	clearnose skate	goosefish	striped bass	striped bass

Most Species Biomass

# Most biomass
biom_table <- top_species %>% 
  arrange(desc(total_bio)) %>% 
  mutate(rank = row_number()) %>% 
  slice_head(n = 3) %>% 
  ungroup() %>% 
  select(-starts_with("total")) %>% 
  pivot_wider(names_from = size_group, values_from = comname) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels)) %>% 
  arrange(survey_area) %>% 
  select(survey_area, rank, S, `S-M`, M, `M-L`, L) %>% 
  group_by(survey_area) %>% 
  gt() %>% 
    tab_header(title = "Most Biomass Overall")

biom_table

Most Biomass Overall
rank	S	S-M	M	M-L	L
GoM
1	acadian redfish	spiny dogfish	spiny dogfish	spiny dogfish	atlantic cod
2	american plaice	white hake	white hake	white hake	white hake
3	silver hake	haddock	atlantic cod	atlantic cod	pollock
GB
1	little skate	spiny dogfish	winter skate	winter skate	winter skate
2	spiny dogfish	winter skate	spiny dogfish	atlantic cod	atlantic cod
3	haddock	haddock	atlantic cod	spiny dogfish	spiny dogfish
SNE
1	little skate	spiny dogfish	spiny dogfish	spiny dogfish	spiny dogfish
2	spiny dogfish	winter skate	winter skate	winter skate	winter skate
3	silver hake	summer flounder	goosefish	smooth dogfish	goosefish
MAB
1	spiny dogfish	spiny dogfish	spiny dogfish	spiny dogfish	spiny dogfish
2	little skate	clearnose skate	smooth dogfish	smooth dogfish	smooth dogfish
3	spotted hake	smooth dogfish	goosefish	goosefish	goosefish

1kg Bin Changes

How has the fraction of total abundance shifted for these weight-bins:

ggplot(even_breaks_summary, aes(est_year, percent_abund, color = season)) +
  #geom_line() +
  geom_point(size = 0.8, alpha = 0.5) +
  geom_ma(size = 1, n = 5, ma_fun = SMA, linetype = 1) +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  facet_grid(size_range*size_group~survey_area, scales = "free") +
  labs(y = "% / Seasonal Abundance", x = "Year")

How is Abundance Distributed?

These figures really should be standardized into some sort of CPUE to account for differences in effort and area.

# Plot the Northeast Shelf changes
even_breaks_summary %>% 
  mutate(size_range = fct_rev(size_range)) %>% 
  filter(survey_area != "Northeast Shelf") %>% 
  mutate(area_title = "Northeast Shelf") %>% 
  ggplot(aes(est_year, total_size_abund, fill = survey_area)) +
  geom_col(position = "stack") +
  scale_fill_gmri() +
  scale_y_continuous(breaks = scales::breaks_pretty(), labels = label_number(big.mark = ",") ) +
  facet_grid(size_range*size_group~area_title*season, scales = "free") +
  labs(y = "Abundance", x = "Year", fill = "Region")

How is Biomass Distributed?

You can really see the jump at 2008 in these figures if you look at Spring, and at the smaller bins for Fall.

Biomass totals below are from biomass-at-length estimates*

# Plot the Northeast Shelf changes
even_breaks_summary %>% 
  mutate(size_range = fct_rev(size_range)) %>% 
  filter(survey_area != "Northeast Shelf") %>% 
  mutate(area_title = "Northeast Shelf") %>% 
  ggplot(aes(est_year, total_size_bio, fill = survey_area)) +
  geom_col(position = "stack") +
  scale_fill_gmri() +
  scale_y_continuous(breaks = scales::breaks_pretty(), labels = label_number(big.mark = ",") ) +
  facet_grid(size_range*size_group~area_title*season, scales = "free") +
  labs(y = "Biomass (kg)", x = "Year", fill = "Region")

Summary of Trends

Let’s start with these pieces, and we can then discuss doing something with the extremes of the size distribution–changes in small vs. large individuals, where we pre-define size bins as you’ve done in slides 7 and 8, or perhaps using an approach similar to Lora’s approach of defining the upper/lower percentiles of size. - KMillz

---

Part 2: Relationships to External Factors

The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.

# # Model Dataset for Wigley Species Subset
wigley_bmspectra_model_df <- read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))

# Use one for plots
wigley_bmspectra_model_df <- wigley_bmspectra_model_df %>% 
  mutate(survey_area = case_when(
    survey_area == "GoM" ~ "Gulf of Maine",
    survey_area == "GB" ~ "Georges Bank",
    survey_area == "SNE" ~ "Southern New England",
    survey_area == "MAB" ~ "Mid-Atlantic Bight"),
  survey_area = factor(survey_area, levels = area_levels_long),
  season = fct_rev(season))

Bottom Temperature Changes

Bottom temperatures are from the DuPontavice and Saba combined bottom temperature product. This data product has been clipped to the study areas and seasonally averaged to match the survey sampling seasons.

# Model linear trends in time
btemp_mods <- lmer(
  formula = bot_temp ~ est_year * survey_area * season + (1|est_year),
  data = wigley_bmspectra_model_df)



# change the stupid labels around because
btemp_preds <- as.data.frame(
    # Model predictions
    ggpredict(
        btemp_mods, 
        terms = ~ est_year * survey_area * season) ) %>% 
    mutate(
        survey_area = factor(group, levels = area_levels_long),
        season = factor(facet, levels = c("Spring", "Fall")))
btemp_signif <- emtrends(
      object = btemp_mods,
      specs =  ~ survey_area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
btemp_trend_predictions <- left_join(btemp_preds, btemp_signif)  %>% 
  mutate(survey_area = group,
         survey_area = factor(survey_area, levels = area_levels_long),
         season = factor(season, levels = c("Spring", "Fall")))




# Plot temperatures
temp1 <- wigley_bmspectra_model_df %>% 
  ggplot() +
  geom_ribbon(
    data = filter(btemp_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    aes(est_year, bot_temp, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, bot_temp, color = season, linetype = "5-Year Smooth"),
    size = 1, n = 5, ma_fun = SMA) +
  geom_line(
    data = filter(btemp_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_color_gmri() + 
  scale_fill_gmri() + 
  #scale_linetype_manual(values = c(1,1)) +
  facet_grid(survey_area~., scales = "free") +
  theme(legend.position = "right") +
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  labs(y = "Bottom Temperature", x = "Year", color = "Season",
       fill = "Season", linetype = "Trend")



temp1 

Bottom temperatures have increased in all four regions (annual averages), but seasonally spring temperatures in SNE and MAB are stable.

Bottom Temperature Distributions

The distribution of seasonal bottom temperatures are such that Gulf of Maine experiences a much smaller difference in seasonal temperature swings than the other regions. In the Gulf of Maine the difference in bottom temperature from Spring to Fall can be smaller than the differences in spring bottom temperature along the whole Northeast Shelf.

# Plot the distribution as a raincloud plot
wigley_bmspectra_model_df %>% 
  ggplot(aes(fct_rev(survey_area), bot_temp, color = season, fill = season)) +
   ggdist::stat_halfeye(
     alpha = 0.4,
     adjust = .5, 
     width = .6, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA
  ) + 
  geom_boxplot(
    fill = "transparent",
    width = .15, 
    outlier.shape = NA
  ) +
  ## add dot plots from {ggdist} package
  ggdist::stat_dots(
    ## orientation to the left
    side = "left", size = 1,
    ## move geom to the left
    justification = 1.12, 
    ## adjust grouping (binning) of observations 
    binwidth = .1) + 
  scale_color_gmri() + 
  scale_fill_gmri() + 
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  coord_flip() +
  theme(legend.position = "none") +
  labs(y = "Bottom Temperature", x = "", color = "Season")

The overlap in bottom temperature distributions is very similar to the overlap in mass spectra exponents. Gulf of Maine spectra and bottom temperatures share a good region of overlap, and the regions that experience larger seasonal fluctuations in BT see larger seasonal fluctuations in their size distributions.

wigley_b_dist_plot <- wigley_bmspectra_model_df %>% 
  ggplot(aes(fct_rev(survey_area), b, color = season, fill = season)) +
  ggdist::stat_halfeye(
     alpha = 0.4,
     adjust = .5, 
     width = .6, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA) + 
  geom_boxplot(
    fill = "transparent",
    width = .15, 
    outlier.shape = NA) +
  ## add dot plots from {ggdist} package
  ggdist::stat_dots(
    ## orientation to the left
    side = "left", size = 1,
    ## move geom to the left
    justification = 1.12, 
    ## adjust grouping (binning) of observations 
    binwidth = .005) + 
  scale_color_gmri() + 
  scale_fill_gmri() + 
  coord_flip() +
  labs(
    y = "ISD Exponent", 
    x = "",
    title = "Individual Body Mass Distribution (1g - Max)")



# Plot next to bottom temperatures
wigley_b_dist_plot

This makes me lean towards belief that the size distribution is more responsive on shorter time scales to the ambient water temperatures in order to align so well with seasonal temperature swings. To me this suggests that as a whole the community is likely better equipped to respond to the long-term trends if they track larger absolute swings within the year.

wigley_bmspectra_model_df %>% 
  ggplot(aes(bot_temp, b)) +
  geom_point(aes(color = season, shape = survey_area)) +
  geom_smooth(method = "gam", color = "black") +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(title = "Broad Relationship to Ambient Bottom Temperature")

Total Landings

The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.

wigley_bmspectra_model_df %>% 
  distinct(est_year, survey_area, total_weight_lb) %>% 
  ggplot(aes(est_year, total_weight_lb/1e6)) +
  geom_col(alpha = 0.5) +
  geom_ma(aes(color = "5-Year Smooth"), size = 1, n = 5, ma_fun = SMA, linetype = 1) +
  scale_color_manual(values = "orange") +
  facet_wrap(~survey_area, ncol = 1) +
  scale_y_continuous(labels = label_number(suffix = "M")) +
  labs(y = "Total Landings (live lb.)", x = "Year", color = "Moving Average")

Zooplankton Community Indices

The ecodata package has changed the zooplankton indices around that it carries. It now contains “regime” metrics for 26 of the main zooplankton taxa.

zoo_regimes <- ecodata::zoo_regime
# zoo_anoms <- ecodata::zoo_abundance_anom
zoo_regimes %>% 
  filter(str_detect(Var, "calfin|chaeto|pseudo")) %>% 
ggplot(aes(Time, Value, color = Var)) +
  geom_point(
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(linetype = "5-Year Smooth"),
    size = 1, n = 5, ma_fun = SMA) +
  facet_grid(EPU~.) +
  labs()

---

Driver Correlations

Driver Modeling

This section is likely to change

The following model has been used to determine the relationship between landings and seasonal bottom temperature on b.

lmer(
  b ~ survey_area*season*scale(roll_temp) + survey_area * scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)

#### Model 2: Temperature & Fishing ####


# Drop NA's
wtb_model_df <- drop_na(wigley_bmspectra_model_df, total_weight_lb, bot_temp, b) %>%
  rename(area = survey_area) %>% 
  left_join(area_df) %>% 
  # Do rolling BT within a season * region
  group_by(survey_area, season) %>%
  mutate(
    roll_temp = zoo::rollapply(bot_temp, 5, mean, na.rm = T, align = "right",  fill = NA)) %>% 
  ungroup() %>% 
  mutate(
    yr_num = as.numeric(est_year),
    yr_fac = factor(est_year),
    survey_area = factor(survey_area, levels = area_levels),
    season = factor(season, levels = c("Spring", "Fall")),
    landings = total_weight_lb,
    yr_seas = str_c(season, est_year))


# Make the model
mod2 <- lmer(
  b ~ survey_area*season*scale(roll_temp) + survey_area*scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)


# summary(mod2)
# check_model(mod2)

Temperature Relationship

# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    mod2, 
    terms = ~ roll_temp*survey_area*season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = mod2,
  ~survey_area * season,
  var = "roll_temp")


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = mod2,
  specs =  ~ survey_area*season,
  var = "roll_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))




# Limit temperature plotting range
actual_range <- wtb_model_df %>% 
  group_by(season, survey_area) %>% 
  summarise(min_temp = min(bot_temp)-2,
            max_temp = max(bot_temp)+2,
            .groups = "drop")



# Plot over observed data
mod2_preds %>% 
  left_join(actual_range) %>% 
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>% 
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = survey_area), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wtb_model_df,
    aes(roll_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Exponent of ISD",
       title = "Wigley Species, Body Mass Spectra",
       x = "5-Year Rolling Average Bottom Temperature")

Landings Relationship

# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    mod2, 
    terms = ~ total_weight_lb * survey_area * season))   %>% 
  mutate(
    survey_area = factor(group, levels = area_levels),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = mod2,
  ~survey_area * season,
  var = "total_weight_lb")
#summary(trend_df, infer= c(T,T))


# Drop effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = mod2,
  specs =  ~ survey_area*season,
  var = "log(total_weight_lb)") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T),
    group = str_c(survey_area, "X", season))




# Plot over observed data
mod2_preds %>% 
  left_join(select(mod2_emtrend, survey_area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wtb_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(survey_area~., scales = "free") +
  scale_color_gmri() +
  scale_x_continuous(transform = "log10", labels = label_log(base = 10)) +
  labs(y = "Body Mass Spectra Slope (b)",
       title = "Wigley Species, Body Mass Spectra",
       x = "Total Landings (lb.)")

Driver Model Summary

# mod2 %>% 
#   gtsummary::tbl_regression()   %>% 
#   add_global_p(keep = T) %>% 
#   bold_p(t = 0.05)  %>%
#   bold_labels() %>%
#   italicize_levels() %>% 
#   as_gt() %>% 
#   tab_header(title = "Wigley Community - Bodymass Distribution & Drivers")

Side Stories

Everything below is digging into different nuances of the approach or the community composition. These things are of secondary importance and will be moved out.

Spiny Dogfish Sensitivity

Spiny dogfish constitute a large majority fraction of biomass in the MAB during the Spring survey. Their numbers have also been increasing over time. This acts to shallow the spectra as they are one of the larger species routinely sampled by the trawl survey.

If we remove them from the estimation and re-run the slopes this is what changes.

How much can one species shift the spectra, and is this a definitive/conclusive proof that its the dogfish shallowing the slope in MAB?

# # Run MAB with and without spiny dogfish
# # Run the year, season, area fits for the filtered data
# no_dogs_bmspectra <- group_binspecies_spectra(
#   ss_input = filter(wigley_in, comname != "spiny dogfish"),
#   grouping_vars = c("est_year", "season", "survey_area"),
#   abundance_vals = "numlen_adj",
#   length_vals = "length_cm",
#   use_weight = TRUE,
#   isd_xmin = 1,
#   global_min = TRUE,
#   isd_xmax = NULL,
#   global_max = FALSE,
#   bin_width = 1,
#   vdiff = 2)

# # Save those
# write_csv(
#   no_dogs_bmspectra,
#   here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv"))

# Read them in, do some plotting
no_dogs_bmspectra <- read_csv(
  here::here("Data/processed/wigley_species_nodogfish_bmspectra.csv")) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels))


# Format a little
no_dogs_bmspectra <- no_dogs_bmspectra %>% 
  mutate(est_year = as.numeric(est_year))


# Join them together
dfish_results <- bind_rows(
   list(
     "Wigley Full" = wigley_bmspectra,
     "Wigley w/o Spiny Dogfish" =  no_dogs_bmspectra), 
   .id = "community") %>% 
  mutate(metric = "bodymass_spectra") %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels))




# Get trends for no dogfish
nodfish_mod <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = no_dogs_bmspectra)
bmspectra_mod_wig <- lmer(
  formula = b ~ est_year * survey_area * season + (1|est_year),
  data = wigley_bmspectra)

# Put predictions
dogfish_sens_predictions <- bind_rows(list(
  "bodymass_spectra-Wigley Full" = get_preds_and_trendsignif(bmspectra_mod_wig),
  "bodymass_spectra-Wigley w/o Spiny Dogfish" = get_preds_and_trendsignif(nodfish_mod)), 
  .id = "model_id") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric))



# Plot the differences
ggplot() +
  geom_ribbon(
    data = filter(dogfish_sens_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = dfish_results,
    aes(est_year, b, color = season),
    size = 0.8, alpha = 0.5) +
  geom_line(
    data = filter(dogfish_sens_predictions, non_zero == TRUE),
    aes(x, predicted, color = season),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(survey_area ~ metric*community, scales = "free") +
  labs(
    x = "Year",
    y = "Exponent of Size Spectra",
    title = "Shift in Spectra Trend(s) with Exclusion of Spiny Dogfish")

Story Thoughts

Below is a table of some papers operating in the same area or around the same topic that are particularly relevant:

shelf_papers <- tribble(
  ~"Author", ~"Year", ~"Conjecture",
  "Friedland et al.", 2024, "Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.",
  "Krumsick", 2020, "Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios"
)


shelf_papers %>% gt()

Author	Year	Conjecture
Friedland et al.	2024	Body size has declined despite more abundance & biomass. Competition is preventing larger sizes.
Krumsick	2020	Empirical slopes from trawl survey in labrador sea were steeper than theoretical from nitrogen stable isotope ratios

There are a couple core ideas about what changes have occurred or have been ocurring over the Northeast Shelf with relevance to the community that is sampled by the trawl survey:

---